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Abstract 

This paper presents an overview of physical ideas and mathematical 
methods for implementing non-smooth and discontinuous substitutions in 
dynamical systems. General purpose of such substitutions is to bring the 
differential equations of motion to the form, which is convenient for fur- 
ther use of analytical and numerical methods of analyses. Three different 
types of nonsmooth transformations are discussed as follows: positional 
coordinate transformation, state variables transformation, and temporal 
transformations. Illustrating examples are provided. 

1 Introduction 

Discontinuities of states in physical models often represent the result of inten- 
tional idealization of abrupt but still smooth changes in dynamical characteris- 
tics. Such idealizations help to skip complicated details of modeling on relatively 
small intervals of the dynamics. From the mathematical standpoint however, 
every discontinuity of states breaks the system and thus increases the system 
dimension as many as twice even though the equations may remain the same 
before and after the discontinuity. This work will outline and illustrate different 
ideas of preventing the dimension increase while dealing with discontinuities of 
system states. Regardless physical principals and mathematical specifics, pre- 
liminary modification of descriptive functions represents a common feature of 
such approaches. Note that many analytical methods dealing with dynamical 
systems preliminarily adapt the equations of motions through different substi- 
tutions and transformations in order to ease further steps of analyses. Although 
a universal recipe for such substitutions is rather difficult to suggest, there are 
some general principles to follow. For instance, classes of transformations should 
comply with classes of systems considered. The term classes remains intention- 
ally unspecified here since it may indicate any generic feature of the model, such 
as linearity or nonlinearity, a class of smoothness, or other mathematical prop- 



erties. In particular, this work focuses on substitutions including nonsmootli or 
discontinuous functions. 

Generally speaking, differential operations with nonsmootli functions require 
generalized interpretations of equalities as integral identities, in other words, in 
terms of distributions |50) . In linear cases, such interpretations are usually quite 
straightforward since distributions represent linear functionals [ST]. Nonlinear 
models however impose certain structural constraints on the presence of non- 
smooth or discontinuous functions in differential equations |14| , |15j . Moreover, 
whether or not some combinations of discontinuous functions are meaningful 
may depend upon physical contents of variables participating in such combi- 
nations |35| . From the mathematical standpoint, physical interpretations allow 
for narrowing families of smooth functions that have discontinuities in their 
asymptotic limits. As a result, some combinations of discontinuous functions 
may acquire certain meanings of distributions. To this end, some introduc- 
tory remarks on transitions to nonsmooth or discontinuous limits in nonlinear 
cases will be introduced. For that reason, consider, for instance, the follow- 
ing product of the Heaviside' unit-step function and the Dirac' delta-function, 
9{t — l)S{t — 1). Generally speaking, such a 'product' is undefined due to the 
discontinuity of 6{t — 1). In other words, no certain number can be prescribed 
to the integral 9{t — l)6{t — l)dt, unless it is known what kind of processes 
both functions actually represented before their transition to discontinuous lim- 
its. For illustrating purposes, let 6e{t - 1) = e-2cosh-2[(t - l)/e'^]/2 be the 
force applied to a unit mass whose velocity is zero at i = — oo. Then the velocity 
is described by v{t) — 9^{t—l) = (l + tanh[(i— l)/£^])/2 due to the relationship 
6e(t — 1) — 6e{t — 1); see Fig. 1 for illustration. 

Based on standard definitions, it can be shown that, in terms of weak limits, 
6e{t — 1) gives 6{t — 1) whereas 9^{t — 1) gives 9{t — 1) as e — > 0. Note that 
the above choice for ^^(t — 1) is not unique - there are many other cases leading 
to the same result. Regardless magnitude of e, however, the work done by this 
force is calculated as follows 



/oo 
9,{t-l)S,{t-l)dt = 
-OO 



9,it~l)K{t~l)dt = l[9eit-i)r\°°^ = l 

Therefore, transition to the limit e — > gives 

^oo 11/"°° 

/ 9{t - l)d{t - l)dt = - = - 5{t-l)dt 

J -oo ^ ^ J -oo 

In terms of distributions, the above expression mean^ 

9{t-l)5{t-\) = \6{t-l) (1) 



^Strictly speaking, a so-called testing function must be added to the integrands as a factor 
in order to completely prove this statement. 
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The left-hand side in ([I]) can be interpreted as a mechanical power generated 
by the impact force S{t— 1). To some extent, the result ([I]) is based on the com- 
mon sense that such a power must exist regardless mathematical ambiguity of 
the left-hand side of identity ([ij . Note that interpretation ([T]) becomes possible 
due to the 'constraint' couphng the two functions, 6g{t — 1) = S^{t — 1), which 
is supposed to hold in the limit case e ^ 0, 

e{t ~ 1) ^ S{t - I) (2) 

were the over dot means Schwartz derivative. 

This standard relationship of the theory of distributions as well as its differ- 
ent variations will be used below. Relationship ([iJ may also serve as a basis for 
other definitions, for instance, 

sgait ~ l)5{t - 1) = (3) 

where sgn(t - 1) = 29{t - 1) - 1. 

While using ([I]) or ([3|, both functions 9{t~l) and S{t~l) must be considered 
together with their 'smooth prehistories' generated by the same family of smooth 
functions 5^{t—l). If the limit functions 0(^—1) and 5{t—l) are based on different 
families of smooth functions then relationships ([T]) or (|3| may not hold; see the 
next section for more details. 

To conclude this, as follows from the above remarks, despite of the universal 
notations, discontinuous functions may still inherit some features of the gen- 
erating families of smooth functions. Therefore, analytical manipulations with 
discontinuous and delta-functions must account for both physical content of the 
problem and mathematical structure of equations. In the most direct way, such 
complications can be avoided by considering models on different time intervals 
and introducing appropriate matching conditions for the corresponding pieces of 
solutions. In many cases, however, the matching times are a priori unknown and 
must also be determined from the same matching conditions. As mentioned at 
the beginning, this work gives an overview of another approaches satisfying the 
matching conditions automatically through specific non-smooth transformations 
of variables. Briefly, the paper is organized as follows: 1) Caratheodori equa- 
tions and discontinuous substitutions |14| for systems under external pulses and 
wave propagation problems, 2) non-smooth coordinate space transformations for 
impact systems with elastic perfectly stiff constraints and possible extensions on 
non-elastic constraints (SU], [ZOIj [ZI]> 3) non-smooth state (phase) space trans- 
formation 12V , [24J , '23| and applications to modeling the impact dynamics with 
an arbitrary coefficient of restitution, 4) non-smooth temporal transformations 
|40j . [41] for periodic and non-periodic motions, analytical and semi-analytical 
tools, unstable periodic and chaotic motions, 5) different combined approaches 

m. m- 
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2 Nonsmooth coordinate transformations 



2.1 Dynamical systems with distributions included as sum- 
mands: linear discontinuous transformations 

Different classes of differential equations, linear and nonlinear, with additively 
involved distributions were considered Filippov |14j. In particular, the methods 
of reducing such equations to Caratheodori equations were described. Briefly, 
such methods eliminate distributions from the differential equations based on 
predictions for local dynamical effects of the distributions. Such effects are 
simply included into the corresponding substitutions so that new equations are 
free of singular terms and thus satisfy the so-called Caratheodory conditions. As 
a result, it is possible to prove existence and investigate different properties of 
solutions. Such a generalization is based on the integral form of the differential 
equation x ~ f{t, x) with a continuous right-hand side 



x{t) = x{tQ)+ J f{s,x(s))ds (4) 



If the function f{t,x) is discontinuous in t, but still continuous in x, then 
the functions satisfying Q can be considered as solutions of the equation 
x = f{t,x). Generally, integration in Q should comply with the concept of 
Lebesgue integral. In this case, the function /(i, x) does not have to be point- 
wise definecQ However, the following Caratheodory conditions must be satis- 
fied: In the domain D of the (t, x)-space, the function f{t,x) is 1) defined and 
continuous in x for almost all t, 2) measurable in t for each x, and 3) the esti- 
mate |/(t, < m{t) holds, where m(t) is a summable function on each finite 
interval, if t is not bounded in D. 

As a simplified illustration, let us consider a single degree-of-freedom system 
whose velocity, v ~ v{t), is described by the differential equation 

v + kv^ ^qSit-ti) (5) 

where k, q and ti are constant parameters, and S is the Dirac'delta function. 

Equation ([s]) describes a unit-mass particle in a nonlinearly viscous media 
with the cubic dissipation law. In this case, the (S-input generates a step-wise 
discontinuity of the response v(t) ai t = ti, nevertheless the nonlinear opera- 
tion in ( |5| r emains meaningful. Moreover, the (5-pulse can be eliminated from 
equation ([5j) as follows. Since the velocity must be bounded, then only the in- 
ertia term can compensate the impact force on the right-hand side of equation 
([5|. In other words, near the impact time t = ti, equation ([S]) is linearized as 
V ~ q6 {t — ti). Comparing this with equation ([2|, gives solution 

V = u + q9{t-ti) (6) 



■^Practically, such an extension almost never contradicts to the physical contents of mod- 
eling; recall that the differential equations of motion are derived from variational principles 
formulated in the integral form. 
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where u is an arbitrary constant. 

Let us assume now that m is a function of time u = u{t) such that sohition 
([6| becomes vahd for the original equation ([5]) on the entire time interval, — oo < 
t < oo. Then, substituting ([6| in ([5|, and taking into account that 6{t — ti) = 
d{t- ti) and 0"^ [t - ti) = 9 {t - ti), gives 

u + k[u^ + {q^ + 3uq^ + 3u'^q)0 {t - h)] = (7) 

In contrast to ([s]), equation Q includes no (5- function and thus admits visu- 
alization of its phase flow on the (i, u)-plane with the related qualitative study. 
Note that substitution ^ generahzes on equation 

C30 

f{v,t)+J2Q^Ht-t^) (8) 
i=l 

where f{v,t) is assumed to have no singularities within some domain of the 
(v,t) -plane, and {qi} and {ti} are constants. 

In this case, the following substitution eliminates all the (5-functions from 
equation ([s]) 

oo 

v{t) = u{t) + q,e it - U) (9) 

i=l 

Further generalization on the vector-form equations is quite obvious. Com- 
plications may occur however when the structure of original equations is changed 
as described in the next section. 



2.2 Distributions as parametric inputs in dynamical sys- 
tems 

This specific case illustrates the major issue of using distributions for modeling 
dynamical systems. Following reference |14| . consider the linear initial value 
problem 

v + k5{t-ti)v ^ (10) 
«(0) = 

where k and ii > are constant. 

Since i; = for t < and t > ti, then 

V ^ vo[l - X9 {t - h)] (11) 

where A is yet unknown constant parameter of discontinuity such that v = vq 
for t < ti and v — ~Xvq for t > ti. 



Equation (10) shows that the jump of the solution depends on the behavior 
of solution itself near the point t = ti. In this case, there is no certain choice 
for A since its magnitude depends upon additional assumptions regarding the 
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model [13]. To clarify the details, let us substitute (11) in (10) with the intent 
to find A. This gives 



(fc - x)S{t - h) - k\e{t ~ ti)S{t -h) = 



(12) 



It was mentioned in Introduction that the combination 9(t — ti)S{t — ti) 
has no certain meaning in the distribution theory. However, the form of (111 



together with equation ( 10 ) dictate 



and therefore, 



9{t - ti)S{t - ti) = adit - ti) 



A 



1 



(13) 



(14) 



where the magnitude of parameter a depends on the model' assumptions. 

For instance, the number a = 1/2 was obtained in introductory example (1). 
In the present case, it is unknown whether the functions 9{t — ti) and S{t — ti) 
are generated by the same family of smooth functions. Nevertheless, physical 
approaches to determining a may be similar to that described in Introduction. 
Namely, let us replace the Dirac function 6{t — ti) in equation ( 10 1 by its smooth 
preimage Se{t — ti) defined in Introduction. As a result, equation (10) takes the 
form of a regular separable equation whose solution is 



V = Vq exp 



— k [ tanh — + tanh — ^ — 
2 \ 



This gives 



vo{l-[l-cxp{-k)]0{t-ti)} as e->0 



(15) 



(16) 



Comparing ( 16 1 to ( 11 ) and taking into account ( 14 ), gives A 
and therefore 

1 1 

a = — — 

1 — exp(— fc) k 



exp(— fc) 
(17) 



Expression (171 shows that definition ( |13[ ) depends upon the model param- 
eter k; see Fig. 2 for illustration. In particular, a — 1/2 as fc — 0, and this 
brings us back to the case considered in Introduction. 

Finally, note that substitutions of type ^ and ([9| as well as its different vari- 
ations are widely used in the literature to describe moving discontinuity waves 
[68] , [35] , [17] . In most such cases though there is no explicit source of disconti- 
nuities that naturally occur as a result of wave shape evolutions predetermined 
by initial perturbations and inner properties of wave model. 



2.3 Nonsmooth positional coordinates 

The idea of nonsmooth coordinates associates with elastic but perfectly stiff 
barriers reflecting moving particles in a mirror-wise manner. Since outcome of 



6 



such reflections is predictable then it can be built into the mechanical model in 
advance through the corresponding nonsmooth coordinate transformations. It 
was shown in references |69j . [70j . and |71j that introducing nonsmooth coordi- 
nates simply eliminates barriers by unfolding the configuration space to include 
the area behind the barrier. As a result, the differential equations of motion are 
derived on the entire time interval with no need in formulation of impact con- 
ditions. For illustrating purposes, consider the following A^-degree-of-freedom 
Lagrangian system 



^ = ^E€-Efc.fe+i-*)' (18) 
hm < 1 (19) 

qo{t) = qM+l{t) = (20) 

This is a chain of unit-mass particles connected by linearly elastic springs of 
stiffness ki . Perfectly stiff elastic constraints are imposed on each of the coordi- 



nates according to (19). Although Lagrangian (18) generates linear differential 



equations, these equations alone do not describe the entire system. Due to the 



presence of constraints (19), the system is actually strongly nonlinear that be- 
comes obvious in adequately chosen coordinates. Transition to such coordinates 
is described by 

q, = r(x,) (21) 
where r is the triangular sine-wave 



2 . . nx 

T (x) — — arcsmsm— - = 

TT 2 

T (x) = r (4 + x) 



X for — 1 < X < 1 

-x + 2 for 1 < a; < 3 



(22) 



Note that notation (22 1 and normalization of the period differ from those in- 
troduced in original works [SS] , [7D] , [7T] . The only reason for such modification 
is to deal with the triangular wave of unit slop^ 



[r'{x)f 



(23) 



for at least almost all x. 

The coordinate transformation (21) brings system (18) through (20) to the 
form 



L 

xo{t) 



1 



N 



N 



i = l 2=0 

XN+l{t) = 



h[T{Xi+l) - T{Xi)Y 



(24) 
(25) 



■^Although this condition does no matter for the method introduced in references [69) . |70) . 
|71l , Section 3 of the present paper describes another method for which normahzation | |23[ l is 
essential. 
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It is seen from ( 24 ) that transformation ( 22 ) preserves the quadratic form of 



kinetic energy while the constraint conditions ( 19 ) are satisfied automatically 
due to the property |r(a;)| < 1. In contrast to (18), Lagrangian (24) com- 



pletely describes the model on the entire time interval < i < oo. However, 
in terms of the new coordinates Xi, the potential energy acquired a non-local 
cell-wise structure so that the corresponding differential equations of motion are 



essentially nonlinear; for example, see (26) below. Now every impact interac 



tion with constraints is interpreted as a transition from one cell to another as 
illustrated below on the two degrees-of-freedom model, N = 2. In this case, 
Lagrangian (24) gives the differential equations of motion on the infinite plane 
— oo < a:^ < oo (i = 1, 2) with no constraints 



xi + [{ko + ki)T{xi) ~ kiT{x2)]T'{xi) = 
X2 + [{kl + k2)T{x2) - kiT{xi)\T'{x2) = 



(26) 



Fig. 3 shows the corresponding equipotential energy levels and a sample 
trajectory of beat-wise dynamics represented by Figs. 4 and 5 in the original 
coordinates. Fig. 3 shows, for instance, that the system is trapped in some 
cells for the energy exchange process. After one of the two masses accumulated 
the energy, which is sufficient to reach the barrier, the impact event happens 
accompanied transition to another cell. The fact of energy exchange inside a 
trapping cell is confirmed by the transversality of incoming and outcoming pieces 
of the trajectory. As long as the mass remains in impact regime, its trajectory 
is passing through one cell to another until the system is trapped again in 
rather another cell for a new energy exchange process. A similar geometrical 
interpretation but for impact mode dynamics was introduced earlier in |65j . 
where the impact modes were associated with 'hidden geometrical symmetries' 
revealed by periodic patterns of equipotential lines as shown in Fig. 6. 

In particular, closed form analytical solutions for different impact modes 
were obtained by means of the averaging procedure. Note that, according to the 
original works [51] , [70] , [7T] , applicability of the averaging procedure constitutes 
the major advantage given by transformation (21) since infinite impact forces 
are effectively eliminated from the system. 

Similar kind of visualization for a two-degree-of-freedom vibrating system 
with only one mass under two-sided constraint condition was used in I46j . 

Let us consider a one-degree-of-freedom unit mass oscillator with the poten- 
tial energy of restoring force P (q), 



(27) 



whose motion is limited by the interval 



-l<q (28) 
It is assumed that the oscillator collides with a one-sided perfectly stiff bar- 



rier ai q — ~1 with no energy loss. In this case, the constraint (28) is eliminated 
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by the space unfolding coordinate transformation g — > a:: 



g = -l + |2; + l| (29) 

where the constant shifts in ( 29 ) are chosen to preserve the origin and barrier 
positions. 



Substituting (29) in (27), gives 

L^]^±^-P{\x + l\-l) (30) 

The system' phase trajectories are described by the energy integral i;^/2 + 
P(|a; + l| — 1) = Const. This family of curves is illustrated by Fig. 7 on the phase 
plane xx at different energy levels in the case P{q) = g^/2. In particular, it is 
seen that high-energy trajectories are non-smooth due to mirror-wise reflections 
against the barrier dX x — ~1 {q — —\). As a result, the phase portrait of the 
oscillator has effectively non-local structure, where any transition from one side 
of the plane to another is the result of impact interaction with the barrier. 

Although original works [53], [ZD]) |ZI]) deal with illustrating models of de- 
terministic dynamics, further applications were shifted mostly into the area of 
random vibrations may be due to earlier works ^9j and 12J. Different analytical 
and numerical tools in this area were developed during recent few years; see 
references [3H], (TU], [H], [5T]. From the standpoint of practical applications, in- 
elastic effects of interactions with stiff constraints become essential. Generally, 
impact dissipation effects can be modeled by the dissipative term [71], [4], [11] 
(1 — k)x\x\5-{x), where 5-{x) is a specific rule rather than conventional Dirac 
function. According to this rule, the impulsive damping acts right before the 
result of such damping namely velocity jump occurs. Such damping model is 
justified if the restitution coefficient k is close to unity so that the factor 1 — k is 
small. In this case, the integral effect of the impulsive damping can play the role 
perturbation within asymptotic procedures, in which the velocity x is given by 
an unperturbed system and therefore remains continuous. Distributed viscous 
damping effects still can be described in a regular way by continuos terms as it 
was done, for instance, in 3J. Non-elastic impact interactions with constraints 
can be modeled also in a purely geometrical way under some conditions on the 
class of motions though [71] . 

Another approach that deals with the class of nonsmooth 'non-conservative' 
transformations is described in the next subsection. 

2.4 Nonsmooth transformation of state variables 

Generalized approaches to eliminating non-elastic constraints should obviously 
involve both types of the state variables - coordinates and velocities. For il- 
lustrating purposes, let us consider the case of harmonic oscillator under the 
constraint condition 

X = Ax (31) 
> (32) 
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where x — \x\(t), xi(t)Y is the system' state vector such that X2 — ii, and 



A = 



1 



(33) 



It is also assumed that every colhsion with the constraint at xi = is 
accompanied by a momentary energy loss characterized by the coefScient of 
restitution k, in other words 



X2(t* +0) = -KX2it*-0) 



(34) 



where t* is the collision time, at which xi{ t*) = 0. 

The idea is to unfold the phase space in such way that the energy loss occurs 
automatically whenever the system crosses preimage of the line xi — 0. The 
corresponding non-conservative transformation was introduced in |21) , |24j , [2 3) 
as a transformation of state vector, x — > y, of the form 



x = Sy 



(35) 



where y = [s{t), v{t)]'^ is a new state vector, and the transition matrix is given 

by 

' 1 

1 — fcsgn(sw) 

where k = (1 — k)/(1 + k). 



sgn(sj 



(36) 



Note that transformation (|35| is strongly nonlinear due to the nonsmooth 

(37) 



dependence S = S(y,fc). Nevertheless, substitution (35 1, gives equation 

y = (S-iAS)y 



There is some 'hidden issue' with substitution (35) since the result (37) has 



the same form as it would have in the case of constant matrix S. In fact, the 
matrix S is constant just almost everywhere. A formal substitution of ( [SS] ) in 
(31) would eventually impose specific conditions on distributions similar to ([3]). 

In the component-wise form, expressions (35) and (37) are written as, re- 
spectively. 



xi = xi{s,v) = ssgii{s) 

X2 — X2{s,v) = sga{s)[l — ksgn{sv)]v 



and 



[1 — ksgn{sv)]v 

-uj^s[l + ksgn{sv)]/{l - k^) 



(38) 



(39) 



Now both unknown components of the state vector are continuous, whereas 



effects of non-elastic collisions ( 34 ) are captured by transformation ( 38 ) 
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Finally, consider the general case of one-degree-of-freedom nonlinear oscilla- 
tor 



X2 = -f{xi,X2,t) (40) 

whose motion is restricted to the positive half plane xi > by a non-elastic 
barrier at xi = of the restitution coefficient k. 



Applying transformation p8| ) to system ( 40 ) , gives 
s = [1 — fcsgn(su)]w 

V = -/(xi(s,D),a;2(s,w),t)sgn(s)[l + /csgn(sw)]/(l-fc2) (4;^) 

Although the technique is illustrated on a one-degree-of-freedom model, sim- 
ilar coordinate transformations apply to multiple degree-of-freedom systems by 
choosing one of the coordinates perpendicular to the constraint. For that rea- 
son, it is convenient to use the Routh descriptive function whose normal to the 
constraint coordinate is Lagrangian whereas other coordinates and associated 
momenta are Hamiltonian [24 1. 

Finally, note that adjustments of classes of smoothness of dynamical sys- 
tems by eliminating infinite discontinuities extend the set of applicable analyti- 
cal tools. For instance, problems of stability and bifurcation analyses of impact 
motions were considered in [53] by means of the linearization technique devel- 
oped for dynamical systems with nonsmooth right-hand sides . In particular, 
the fundamental matrix of explicit form and the corresponding characteristic 
equations were obtained. As a result, investigation of stability and bifurcations 
became possible by using conventional methods. In fact, before the transfor- 
mation, discontinuities of phase trajectories as those shown in Fig. 9 would 
complicate any local analyses. However, the transformation improves the class 
of smoothness to the extent which is needed to build major objects of local 
analyses and averaging tools. A regular approach to stability and bifurcation 
analysis in impact systems was proposed in [22j. In particular, it was shown 
that the discontinuous bifurcation of grazing impact can be regularized. This 
may lead to a new interpretation of grazing bifurcations. Namely, after such 
bifurcation, some periodic motion might survive and even preserve stability. 



3 Nonsmooth temporal arguments 

In this section, we describe nonsmooth substitutions of the independent variable, 
which is temporal argument in the present text. It will be shown below that such 
nonsmooth substitutions associate with some common temporal symmetries of 
motions regardless of types of systems. This approach was originally developed 
for the class of strongly non- linear but smooth oscillators [30], [JT]. However, 
its main specifics are clearly seen even in the case of non-oscillatory motion of 
a classic unit-mass particle under returning potential force as shown in Figs. 
10 and 11. The idea is to employ most elementary macrodynamical processes 
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whose specifics nevertheless provide sufficient conditions for observing strong 
nonlinearities |48| . Examples of such 'elementary strongly nonlinear processes' 
are found among the rigid-body motions as those shown in Fig. 12. However, 
the key question is how to bridge the gap between the classes of smooth and 
nonsmooth motions within the same mathematical formalism. It will be shown 
below on simple examples that nonsmooth substitutions of temporal argument 
may play the role of such a bridge. 



3.1 Positive time 



Let the potential energy P{x) be a smooth function of the time dependent 
coordinate x as qualitatively shown in Fig. 10. Then the differential equation 
of motion of a unit-mass particle under the force f{x) = P'{x) is given by 



+ fix) = 



(42) 



The initial conditions are x — xq > and x — vq < at t = to < as shown 
in Figs. 10 and 11. 

As the particle reaches a turning point at some time t = a it makes a U- 
turn. Since equation (42 1 admits the group t — > —t, the reverse motion will 
be symmetric with respect to the U-turn point t = a. Such a 'prediction' builds 



into the differential equation of motion ( 42 ) through the new temporal argument 
t — > s: 

s — \t ~ a\, X — xi^s) (43) 

The mechanical model generating such time argument by its natural motion 
is shown in Fig. 12 (a). Since at almost all i, 



(44) 



then, substituting (43) in (42), gives 



ds'-' 



+ /(.t) = 0, s>0 



under condition 



dx 
ds 



if 



(45) 



(46) 



It is easy to see that the boundary condition (46) eliminates singularity 



of second time derivative that formally occurs due to the non-smoothness of 



substitution (43) 



As follows from Fig. 11, substitution (43) reverses the time direction exactly 



when the particle makes a U-turn. Although the form of the equation remains 
the same, substitution (43 1 gives certain advantages from both physical and 



mathematical standpoints. First, equation [45] even after being drastically sim- 
plified as (P'xjds^ ~ 0, still preserves the major dynamical event, which is the 
U-turn of the particle. In this degenerated case, the smooth potential barrier is 
effectively replaced by a perfectly stiff one as follows from the general solution. 
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X = As{t) + B, where A and B are arbitrary constants of integration. Now, if 
some perturbation series converges for s > 0, then it is automatically converges 
for the entire interval of the original time, — oo < t < oo. 

As another example, let us consider the case of impulsively loaded single 
degree-of-freedom system, 

X + f{x) = ps (47) 
where s — 2d{t — a) and p = const. 



Substituting ( 43 ) in ( 47 ) and taking into account ( 44 ) , gives 



fx 



+ fix)= [P- 



dx 
ds 



(48) 



Eliminating the singularity ps in ( 48 1 , gives the same equation ( 45 1 however 



under non-homogeneous boundary condition 

dx 



ds 



— p if s — 



(49) 



Since substitution (43) is non-invertible in the form t — t{s) on the entire 
time interval, then using the argument s in general case of dynamical system ap- 
pears to be less straightforward but nevertheless possible based on the following 
identity 

t = a + ss (50) 



Due to relationship (44), the combination (50) represents a specific complex 



number with the basis {l,s} |45| . In contrast to the conventional elliptic com- 



plex algebra, the operation 1/t with (50 1 may not hold. Interestingly enough, 
such algebraic structures has been known for quite a long time [8] , [61] mostly as 
abstract mathematical objects with no relation to nonsmooth functions or im- 
pact systems. In the modern mathematical literature, these numbers are often 
referred to as a simple example of Clifford algebras under the name 'hyperbolic 
algebra' |56j with very many synonyms though [1 . Some areas of physics are 
linked to this algebra quite closely [19] , but constructive applications are rather 
limited. 

Let us formulate useful algebraic properties of hyperbolic numbers adapted 



to the form (50 1. For instance, the hyperbolic conjugate to (50) is introduced 
as = a — ss, and the hyperbolic modulus is calculated as \t ^ = \l\tt~ \ = 
\/\o? — s'^\. The term hyperbolic associates with the fact that the relationship 
\t\}^ — p with some fixed p describes a four branched hyperbola on the hyperbolic 
plane with the basis {1, s}: 



±p{cosh( 
±p(sinh C) 



' s sinh ( 
s cosh ( 



±pexp(s(j 
±pexp(s(j 



arctanh(s/a) 
arctanh(a/s) 



The hyperbolic numbers create isomorphism with symmetric 2 x 2-matrixes, 
so that, for instance 



(a + ssf 



a s 
s a 



(51) 
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Note that, in our case, the hyperbohc structure is generated by time depen- 
dent transformations rather than imposed on abstract elements by mathematical 
definitions. As a result, differential and integral properties of hyperbolic num- 
bers can be introduced |45j . So, for practically any function x{t), it can be 
shown that 

x{t)^x{a + ss) ^X{s)+Y (s) s (52) 

where 



X (s) = -[x {a + s) + X {a — s) 
Y (s) = - [x (a + s) — X [a — s) 



Then, taking into account (44), gives 



x{t) = ¥' {s)+X' {s)s+ps 



(53) 



(54) 



where 



Y{0)=p 



(55) 



If x{t) is continuous then obviously p — and derivative (54) also belongs 



to the hyperbolic algebra, otherwise the quantity p may serve for elimination of 



singularities as seen from equation ( 48 ) 



As a simple application, consider the initial value problem 

X + Xx ~ 2pS{t — a) = ps 
x{0) = 



(56) 



where A and p are constant. 



By considering functions X (s) and Y (s) in (52) as new unknowns, then 



pendence of basis {1, s}, gives 



substituting ( 52 ) through ( 54 ) in ( 56 ) , and taking into account the linear inde- 



Y' + XX = 
X' + XY = 
X{a)~Y{a) = 



(57) 



under condition ( 55 ) 



The boundary value problem ( 55 ) is free of (5-functions and admits an obvious 
solution , X = Y ~ pexp(— As). Substituting these X- and K-components in 



( 52 ) , gives solution of the initial value problem ( 56 1 in the form 



X — pexp(— As)(l + s) 



(58) 



The link established between nonsmooth temporal substitution ( 50 ) and the 



hyperbolic structure essentially facilitates analytical manipulations with differ- 
ential equations. For example, instead of using the standard basis {l,s} it is 
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possible to introduce the so-called idempotent basis, associated with the two 
isotropic lines which separate the hyperbolic quadrants as follows 



;(1±^) 



(59) 



= i±, i^i^ = 
Transition to the idempotent basis in ( 50 ) gives 

t = (a + s)«_|_ + (a — s)i_ 



(60) 



Due to the properties of basis (591, substitution (60) possesses 'functional 
linearity.' In other words, 



e = {a + s)°'i+ + {a~ s)°'i^ 

for any real a, and generally, 

x{t) — x[(a + s)z+ + (a — 

= x{a + s)«_|_ + x{a — s)i_ 
= X+{s)i+ + X^{s)i^ 



(61) 



(62) 



The advantage of the idempotent basis is that the differential equations of 
motion with respect to components and X_ are decoupled. Instead the 
boundary conditions become coupled though. Different examples of solutions 
in the idempotent basis are given in recent publications |49j and |48| . 



3.2 Triangular sine-wave time substitution 



Since any vibrating process is a sequence of U-turns then the corresponding 



nonsmooth time substitution can be combined of functions given by (431 with 
different signs and temporal shifts. In periodic case of the period T — A, such 



combination is given by the sawtooth function ( 22 ) , whose argument is replaced 
by time, r — T{t). A mechanical model generating such time substitution by its 



free motion is shown in Fig. 12 (b) while the periodic version of identity (50) is 



t = l + (r-l)- 



if - 1 < t < 3 



(63) 



where = 1, and therefore (63) is a hyperbolic number with the basis {1, t}. 
In physical terms, it follows from (63) that any periodic process, whose 



period is normalized to T = 4, is uniquely expressed through the dynamic 
states of standard impact oscillator in the form [35] 



x{t)=X (r) + r (t) T 



(64) 



where 



X{r) 
Y{t) 



[x (t) +x{2^ t)] 
[x{t)^x{2-t)] 



(65) 
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Identity (64 1 means that the triangular sine and rectangular cosine waves 
capture temporal symmetries of periodic processes regardless specifics of indi- 
vidual vibrating systems. 



Introducing a slow temporal scale in (64) extends the area of applications 
on modulated vibrating processes. In such cases, two variable expansions [2 7) 
can be used by considering the triangular sine time as a fast scale |3H] ■ 

Different applications of nonsmooth argument substitutions with the related 
techniques to problems of theoretical and applied mechanics can be found in [47] , 
[30], [32], [3l], [16], [63], [29], [28], [55], [57], [53], [54], [52], [59], [37], [36], [60], 
[2D] . The methodology was adapted also to the nonlinear normal mode analyses 
and included in monograph [M]. While the idea of NNMs is effective in case 
of weak or no energy exchange, the concept of the limiting phase trajectories 
[30] considers the opposite situation namely intense energy exchanges between 
weakly coupled oscillators or modes [32], [34], [33]. In this case, nonsmooth 
time substitutions are invoked by the temporal behavior of phase angle, which 
is responsible for energy distribution. This resembles the sawtooth wave as the 
energy swing reaches its asymptotic limit. 

Very recently, a class of strongly nonlinear traveling waves and localized 
modes in one-dimensional homogeneous granular chains with no precompres- 
sion were considered in \5S\ . In particular, the asymmetric version of identity 



(64) [42 was applied. As a result, the authors developed a systematic semi an- 
alytical approaches for computing different families of nonlinear traveling waves 
parametrized by spatial periodicity (wave number) and energy. 

3.3 Quasi linear asymptotics for nonsmooth perturbations 

There are many physical and mechanical models described by linear oscillators 
with small but nonsmooth perturbations. The related examples of mechanical 
oscillators were considered in reference [31]. Nonsmooth but still continuous 
characteristics of elastic forces often occur when modeling the dynamics of elastic 
structures with cracks [7], [1], [HS]- More references and mathematical remarks 
on different cases of nonsmooth perturbations can be found in reference . As a 
non-conservative case of nonsmooth perturbations, let us mention the following 
modification of Van der Pol's oscillator [^S] , [TS] 



q + ei\q\-l)q + q = (66) 



where e is a small parameter. 



Note that, using the nonsmooth term in (66), reduces the degree of nonlin- 



earity as compared to the classical Van der Pol's model. A periodic limit cycle 



solution of equation (66 1 was obtained in [18: by taking into account the identity 
\q\ —sgn(q)q in order to facilitate trigonometric expansions for the method of 
multiple scales. The following generalized model was considered in j5j 

q + e{\q\-l)q+{l+ea)q:^eXsint (67) 

where tr is a detuning parameter. 



16 



Recall that classical methods of asymptotic integration require perturbations 
to be as smooth as needed for deriving solutions of certain asymptotic order. 
In such cases, non-smooth argument substitutions may facilitate formulations 
for high-order asymptotic approximations. Moreover, firts-order asymptotic 
solution are obtained exactly in a closed form. Let us illustrate this remark on 
the one-degree-of-freedom oscillator 



= eoj'^9{q)q 



(68) 



where 0{q) is the Heaviside unit-step function. 



The perturbation on the right-hand side of oscillator (68) is a continuous 



but non-smooth function of the coordinate q. Nevertheless, let us represent the 
first-order asymptotic solution in the form 



q = A cos ip + eqiif) + 0{e'^) 

ip = W(l+£7i +0(£2))t 



(69) 



where A = const. > 0, and qi and are yet unknown corrections to the 
generating solution obtained under the condition e = 0. 



Substituting ( 69 1 in ( 68 ) and matching first-order terms of e on both sides 



of the equation, gives 



+ qi = Acosp[2^^ + 0{cos'p)] 



(70) 



where the amplitude A in the argument of unit-step function 9 has been ignored 
as a positive factor of no influence on the output. 

For further comparison reason, let us reproduce first solution of equation 
(70) in terms of trigonometric expansions. According to the idea of Poincare- 



Lindstedt method, the parameter is determined by using Fourier series with 
respect to tp and eliminating the resonance term on the right-hand side of equa- 
tion (70) that gives = —1/4. The first-order approximation qi{p) is obtained 



then in the form 



q 



= A 



cos (p 



l-|+0(e^) 



■ cos 2ip - 



225 



cos Aip 



+ Oie') (71) 



The error of solution ( 71 1 depends upon both the magnitude of e and the 



number of terms retained in the Fourier series. 



Now let us solve equation (70) by applying identity (64 1 to the right-hand 



side of equation ( 70 ) and its unknown 27r-periodic solution. In other words, let 



us represent solution in the form 



qi = X(r)+r(r)e 
T — t{2lp/'k) 

e = e{2ip/Ti) = dT(2ip/Ti)/d{2Lp/n) 



(72) 
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where X and Y are new unknown functions of the argument r, and obviously 
= 1 for almost all cp. 



Then, substituting i72l in ( 70 ) and taking into account the identities cos if 



cos(7rr/2)e and 9{cosip) = (1 + e)/2, gives 



2V d^X 



2V d^Y 



1 TTT , 



Y 



1 \ TTT 

- + 27ijAcos— , Y\r=±l 



(73) 

(74) 



Recall that boundary conditions in (73 1- (74 1 eliminate the singularities of 
differentiation caused by the formal presence of non-smooth functions in ( 72 1 . It 



can be verified by inspection that problem ( 73 ) always admits a solution whereas 
problem (74) is solvable only under the condition ji = —1/4. Note that this 
number is identical to that eliminates the resonance term on the right-hand side 
of equation (70). In the present procedure, however, the number "fi = —1/4 



eliminates the 'imaginary' component from representation (72 ) due to the trivial 



solution of boundary value problem (74) whereas the 'real' component is easily 
determined from (73). As a result, first-order asymptotic solution (69) takes 
the closed form 



A 



cos (f - 



2 cos ■ 



TTTX 



(75) 



Note that solution ( 75 ) appears to have the so-called 'secular term' r sin(7rT/2), 
which is bounded and periodic however with respect to the original temporal 
argument t. 

In order to compare solutions ( 71 ) and ( 75 ) , let us fix the amplitude pa- 
rameter as A = 1.0 in ( 75 ) and then determine the corresponding amplitude 
parameter of solution (71 ) in order to achieve the same initial value q(0). Then, 
selecting e = 0.25 and keeping the three terms of Fourier expansion of order e 



as shown in (71), gives a very good match of solutions (71) and (75) in terms 



of the coordinate q. This is due to quite a rapid convergence of the Fourier 
series in (71) as seen from its coefficients. However, since differentiation slows 



down the convergence then some mismatch in accelerations occurs between the 
two solutions as shown in Fig. 13. In the diagram, the dashed line represents 
numerical solution produced by the built-in Mathcmatica® solver. 



3.4 Nonsmooth time decomposition 

In this subsection, the original temporal argument t is replaced by a sequence 
of nonsmooth arguments {si} {i = 1,2,...) running within the same standard 
interval < < 1. This brings advantage of using bounded time arguments 
and, as shown below, provides a convenient description of impulsively loaded 
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dynamical systems. As a mathematical basis, consider the ramp function, 

sit;d)^^-id+\t\-\t-d\) (76) 

and its first order generahzed derivative, s(t;d), with respect to the temporal 
argument, t; see Figs. 14 and 15, respectively. 

Such kind of functions are quite common for signal analyses p5J. In our 
case, however, physical interpretation of these functions represented in Fig. 12 
(c) are employed. Namely, the function s {t, d) describes positions of a very 
small perfectly stiff bead located initially at the origin x = as shown by the 
dashed circle. The initial time instance t = associates with the event when 
this bead is struck on the left by the identical bead with the velocity v = 1. 
Due to the linear momenta exchange, the reference bead starts moving by the 
law X = x{t,d) until it hits the third bead located a.t x = d = 1. Now, let us 
consider an infinite chain of perfectly stiff identical beads on a straight line at 
the points Xi (i = 0, 1, ...). As soon as no energy loss is assumed, any currently 
moving bead has the same velocity, v — 1. As a result, the linear momentum is 
translated with the constant speed, whereas every bead moves only in between 
its two impact interactions with the neighbors. The motion of the ith bead is 
described by 

(t) = sit- U, di) = ^ {d, + \t-U\-\t~ U+i\) (77) 

where ti is the first interaction time of the ith bead, and di — U^i — U is the 
time interval between its two consecutive interactions with its neighbors. 

Due to the unit velocity, the variable Si can play the role of 'local' time 
associated with the bead moving within the interval Xi < x < Xi^i while the 
'global' time t runs together with the linear momentum across the entire chain 
of beads. However, for any sequence of time instances, A = {tQ,ti, ...}, the 
global time, t £ [to,oo), can be expressed through the sequence of local times, 
{si}, in the form [ii] 

oo 

< = ^ (i» + s,) h (78) 

i=0 

where the basis {si} obeys the following table of products 

SiSj SiS^j ('^^) 

As a result, the following 'functional linearity' property holds for quite a 
general process x{t) 

Coo \ oo 

{U + Si)si\ = ^ a; {U + Si) Si (80) 
1=0 / i=0 

As mentioned at the beginning of this subsection, the area of applications 



of time substitution ( 78 ) includes impulsively loaded dynamical systems. In 
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particular, taking into account (801 and differentiation rules [U], eliminates the 
external impulses from dynamical system, 



X = i{x,t) + J2p^Ht-U), x(i)ei?" (81) 



1=0 



X EE 0, t <to 

where f (x, t) is regular vector- function and Pi are vectors characterizing mag- 
nitudes and directions of the impulses. 

The new system is considered then in the local (bounded) time arguments 
{si} by using general analytical tools such as Lie series expansions. Such ap- 
proach was applied to the Duffing oscillator with no linear stiffness under sine 
modulated random impulses |44) . 

oo 

i + (x + x^ ^ B suites {t~U) (82) 

i=Q 

where C is a constant linear damping coefficient, and B is the amplitude of 
modulation. 

The distance between any two sequential impulse times is given by 

where 77,^ are random real numbers homogeneously distributed on the interval 
[— 1, 1], and /3 is a small positive number, < /3 <C 1. 



Introducing the state vector x = {x,x)'^ — {xi,X2)'^, brings system (82) to 



the standard form (81), where 



f(x) 



X2 \ ^[0 



Note that oscillator (82) is a modification of the well known oscillator, x + 
(x + x^ = Bsint, considered by Ueda [62] as a model of nonlinear inductor in 
electrical circuits. In particular, the result of work |62j . as well as many further 
investigations of similar models, revealed the existence of stochastic attractors 
often illustrated by the Poincare diagrams. Similar diagrams obtained under 
irregular snapshots can be qualified as 'stroboscopic' diagrams. The results of 
the computer simulations [44j show that some irregularity of the pulse times 
can be used for the purposes of a more clear observation of the system orbits in 
the stroboscopic diagrams subjected to a chain of significant structural changes 
as the parameter B increases. When repeatedly executing the numerical code, 
under the same input conditions, such a small disorder in the input results some 
times in a less noisy and more organized stroboscopic diagrams. However, such 
phenomenon itself was found to be a random event whose 'appearance' depends 
on the level of pulse randomization as well as the number of iterations. 
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4 Concluding remarks 



This work outlined very different ways to modeling dynamical systems with 
discontinuities by choosing proper spatial coordinates or temporal arguments 
within the class of nonsmooth functions. Note that nonsmooth transformation 
of state vector ( 38 1 can be viewed as a generalization of nonsmooth coordinate 
transformation (29). Whereas transformation (29 1 is conservative, transforma- 
tion ( 38 ) possesses built in energy sinks pumping out some energy from the 
system whenever it crosses the image of the barrier in the unfolded space. In 
contrast to substitution ^ , both substitutions ( 38 ) and ( 29 1 are essentially 
nonlinear and thus produce essentially nonlinear differential equations of mo- 
tion even when out of constraints the equations of motion are linear. None of 
the above transformation affects the time variable. In contrast, the nonsmooth 
time substitutions described in Section 3 incorporate general temporal symme- 
tries of dynamic processes. Although such symmetries develop most explicitly 
through natural motions of elementary impact models as shown in Fig. 12, the 
corresponding time substitutions impose no constrains on class of smoothness 
of those systems to which such substitutions are applied. One specific feature of 
the nonsmooth time substitutions is that they induce the hyperbolic structure 
of spatial coordinates. Further details on the related mathematical properties 
and physical interpretations are described in 48^ . Finally note that it is possible 
to combine different transformations described in the present survey whenever 
technical reasons for such combinations are present; some examples can be found 
in ESI and EHl- 
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FIGURE CAPTIONS 

Fig. 1: Temporal shapes of the force apphed to the unit mass 6^ and the 
corresponding velocity for £ = 0.5. 

Fig. 2: Impulse produced by the 'product' of unit-step and Dirac delta 
functions versus the model parameter. 

Fig.3: Equipotential energy levels in the unfolded configuration plane and 
a sample dynamic trajectory obtained under the initial conditions at f = 0: 
x\ = 0.5, X2 = 0.0, x\ = 1.0, and X2 = 0.0. 

Fig. 4: The beat-wise impact dynamics of the first mass in its original 
coordinate. 

Fig. 5: The beat-wise impact dynamics of the second mass in its original 
coordinate; the energy exchange between the two masses is seen from the phase 
shift in their time histories. 

Fig. 6: The impact mode trajectories in the unfolded configuration plane: 
in-phase and out-of-phase modes (the diagonal lines), and local modes (the 
horizontal and vertical lines.) 

Fig. 7: A family of phase trajectories of the oscillator with one-sided barrier 
at different energy levels for the case P{q) = 9^/2. 

Fig. 8: The phase trajectory of inelastic impact oscillator in the auxiliary 
coordinates. 

Fig. 9: The original phase plane of the harmonic oscillator with a perfectly 
stiff but inelastic one-sided barrier. 

Fig. 10: A classic particle reflected by smooth potential barrier. 

Fig. 11: The nonsmooth temporal argument s and the particle' motion, 

x{t). 

Fig. 12: Three basic impact models generating nonsmooth time substitu- 
tions: (a) positive time, (b) triangular sine-wave time, and (c) nonsmooth time 
decomposition. 

Fig. 13: Acceleration curves: Poincare-Lindstedt method (thin line), nons- 
mooth time substitution (thick line), and numerical (dashed line); w = 1.0, and 
£ = 0.25. 

Fig. 14: The unit slope ramp function at d = 1.0. 
Fig. 15: First derivative of the ramp function 



28 



3 



-1 



Figure "1: Temporal shapes of the force applied to the unit mass ig and the 
corresponding velocit}' 9^ for e = 0,5, 




Figure 2; Impulse produced by the 'product' of unit-step and Dirac delta func- 
tions versus the model parameter. 
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Figure 3; Equipotential energy le\-el3 in tlie untolded configuration plane and 
a sample djaiamic trajectory obtained under tlie initial conditions at t = 0: 
xi = 0.5, X2 = 0.0, ±1 = 1.0, and ±2 = 0.0. 
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Figure 4: The beat-wise impact dynamics of the first mass in its original coor- 
dinate. 
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Figure 5: The beat- wise impact dynamics of the second mass in its original 
coordinate: the energy exchange between the two masses is seen from the phase 
shift in their time histories. 
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Xi 

Figure 6: The impact mode trajectories in the unfolded configuration plane: 
in-phase and out-of-phase modes (the diagonal lines), and local modes (the 
horizontal and ^'el•tica] lines.) 
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Figi.u e 7; A family of pliase trajectoi'ies of the oscillator with one-sided barrier 

at different energy levels foi' the case P(qj = 
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Figui e S: The phase trajectory of inelastic impact oscillator in the auxiliary 
coordinates. 
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Xi 

Figure 9: The original phase plane of the harmonic oscillator with a perfectly 
stiff but inelastic one-sided barrier. 
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Figure 11; The nonsmooth temporal argument s and the particle' motion, x[i}. 
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Figure 12: Three basic impact models generating nonsmooth time substitu- 
tions: (aj positive time, (hj triangulai sine-wave time, and (c) nonsmooth time 
decomp osition . 
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Figure 13: Acceleration curves: Poincare-Lindstedt method (thin line), non- 
smooth time substitution (tliick line}, and numerical (dashed line); ui = 1.0, 
and e = 0.25. 
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Figure 14: The unit slope ramp function at d = 1.0 
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Figure 15; First derivative of the ramp function 
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